home *** CD-ROM | disk | FTP | other *** search
/ CD Concept 6 / CD Concept 06.iso / mac / UTILITAIRE / RLaB / toolbox / roots.r < prev    next >
Text File  |  1994-05-11  |  1KB  |  57 lines
  1. //-------------------------------------------------------------------//
  2. //  roots
  3.  
  4. //  Syntax: roots(P)
  5.  
  6. //  Description:
  7.  
  8. //  Find roots of polynomial P of degree n.
  9. //
  10. //  P is a vector (row or column) with n+1 components respresenting
  11. //  polyniomial coefficients in decending powers.
  12. // 
  13. //  Return a column vetor containing n roots of P.
  14.  
  15. //  Example:
  16. //  the following script finds the roots of 
  17. //  x^3 + 3*x^2 + 5 = 0
  18. // 
  19. //  > p=[1,3,0,5]
  20. //   p =
  21. //          1          3          0          5
  22. //  > r=roots(p)
  23. //   r =
  24. //                -3.43 + 0i
  25. //             0.213 - 1.19i
  26. //             0.213 + 1.19i
  27.  
  28. //  See Also: poly
  29.  
  30. //  Tzong-Shuoh Yang (tsyang@ce.berkeley.edu)  5/7/94
  31. //
  32. //-------------------------------------------------------------------//
  33. roots = function(P)
  34. {
  35.    local(c,p,r,z,nz);
  36.  
  37.    if (P.n <= 1) {
  38.       error("Argument of roots() is a scalar, not a polynomial.");
  39.    else if (min(size(P)) != 1) {
  40.       error("Argument of roots() must be a column or a row vector.");
  41.    }}
  42.    p = P[:]';
  43.    // trim leading and trailing zeros
  44.    z  = find(p != 0);
  45.    nz = p.n - z[z.n];
  46.    p  = p[z[1]:z[z.n]];
  47.    // find companion matrix
  48.    c = compan(p);
  49.    // find roots and stuff with trailing zeros
  50.    r = [zeros(nz,1);eign(c).val'];  
  51.    // trim complex part if all real
  52.    if (all(imag(r) == 0)) {
  53.       r = real(r);
  54.    }
  55.    return r;
  56. };
  57.